REPORT  DOCUMENTATION  PAGE 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  inst 
data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  a: 
this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  xrt  Reports  (0704-01 
4302  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  tor  failing 
valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. _ 


AFRL-  SR-BL-TR-00- 

6^ 


1.  REPORT  DATE  (DD-MM-YYYY) 

7/25/00 


2.  REPORT  TYPE 

Final  Technical  Report 


3.  DATES  COVERED  (From  -  To) 
4/1/95  -  6/30/98 


4.  TITLE  AND  SUBTITLE 

Computational  Algorithms  for  Aerodynamic  Analysis  and  Design 

Multidisciplinary  Optimization  -  (Subcontract  title  - 
Aerodynamic  Shape  Optimization  Techniques  Based  on  Control  Theory) _ _ 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 

F4  9620-95-1-0259 


5c.  PROGRAM  ELEMENT  NUMBER 


p  AUTHOR(S) 

Luigi  Martinelli  is  the  PI  for  the  subcontract  that  we  are  submitting. 


5d.  PROJECT  NUMBER 

2304/CS 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Mechanical  and  Aerospace 
Engineering  Department 
Princeton  University 
Princeton,  NJ  08544 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

AFOSR/NA 

801  North  Randolph  Road 
Room  732 

Arlington,  VA  22203-1977 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 
AF0SR 


11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 


1 2.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 

13.  SUPPLEMENTARY  NOTES 


20010102  031 


14.  ABSTRACT 

This  document  serves  as  a  final  technical  report  for  the  AFOSR  award  F49620-95-1-0259 .  It 
reviews  the  formulation  and  application  of  optimization  techniques  based  on  control  theory 
for  aerodynamic  shape  design  in  viscous  compressible  flow.  The  theory  is  applied  to  a  system 
defined  by  the  partial  differential  equations  of  the  flow,  with  the  boundary  shape  acting  as 
the  control.  The  Frechet  derivative  of  the  cost  function  is  determined  via  the  solution  of 
an  adjoint  partial  differential  equation,  and  the  boundary  shape  is  then  modified  in  a 
direction  of  descent.  This  process  is  repeated  until  an  optimum  solution  is  approached. 

Each  design  cycle  requires  the  numerical  solution  of  both  the  flow  and  the  adjoint  equations, 
leading  to  a  computational  cost  roughly  equal  to  the  cost  of  two  flow  solutions. 
Representative  results  are  presented  for  viscous  optimization  of  transonic  wing-body 
combinations . 


15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 
OF  PAGES 

a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

32 

19a.  NAME  OF  RESPONSIBLE  PERSON 

LUIGI  MARTINELLI _ 


19b.  TELEPHONE  NUMBER  (include  area 
code>  609-258-6652 


DTIC 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


Aerodynamic  Shape  Optimization  Techniques 
Based  On  Control  Theory 


Antony  Jameson1  and  Luigi  Martinelli2 

1  Department  of  Aeronautics  &:  Astronautics 
Stanford  University 

Stamford,  California  94305  USA 

2  Department  cf  Mechanical  &  Aerospace  Engineering 
Princeton  University 

Princeton.  N3  08544,  USA 


Abstract 

This  document  is  serves  as  a  final  technical  report  for  the  AFOSR  award  F49620-95- 1-0259.  It  reviews 
the  formulation  and  application  of  optimization  techniques  based  on  control  theory  for  aerodynamic  shape 
design  in  viscous  compressible  flow.  The  theory  is  applied  to  a  system  defined  by  the  partial  differential 
equations  of  the  flow,  with  the  boundary  shape  acting  as  the  control.  The  Frechet  derivative  of  the  cost 
function  is  determined  via  the  solution  of  an  adjoint  partial  differential  equation,  and  the  boundary  shape 
is  then  modified  in  a  direction  of  descent.  This  process  is  repeated  until  an  optimum  solution  is  approached. 
Each  design  cycle  requires  the  numerical  solution  of  both  the  flow  and  the  adjoint  equations,  leading  to  a 
computational  cost  roughly  equal  to  the  cost  of  two  flow  solutions.  Representative  results  are  presented  for 
viscous  optimization  of  transonic  wing-body  combinations. 


1  Introduction:  Aerodynamic  Design 

The  definition  of  the  aerodynamic  shapes  of  modern  aircraft  relies  heavily  on  computational  simulation 
to  enable  the  rapid  evaluation  of  many  alternative  designs.  Wind  tunnel  testing  is  then  used  to  confirm 
the  performance  of  designs  that  have  been  identified  by  simulation  as  promising  to  meet  the  performance 
goals.  In  the  case  of  wing  design  and  propulsion  system  integration,  several  complete  cycles  of  computational 
analysis  folic:  ed  by  testing  of  a  preferred  design  may  be  used  in  the  evolution  of  the  final  configuration. 
Wind  tunnel  *^sting  also  plays  a  crucial  role  in  the  development  of  the  detailed  loads  needed  to  complete 
the  structural  design,  and  in  gathering  data  throughout  the  flight  envelope  for  the  design  and  verification 
of  the  stability  and  control  system.  The  use  of  computational  simulation  to  scan  many  alternative  designs 
has  proved  extremely  valuable  in  practice,  but  it  still  suffers  the  limitation  that  it  does  not  guarantee  the 
identification  of  the  best  possible  design.  Generally  one  has  to  accept  the  best  so  far  by  a  given  cutoff  date  in 
the  program  schedule.  To  ensure  the  realization  of  the  true  best  design,  the  ultimate  goal  of  computational 
simulation  methods  should  not  just  be  the  analysis  of  prescribed  shapes,  but  the  automatic  determination 
of  the  true  optimum  shape  for  the  intended  application. 

This  is  the  underlying  motivation  for  the  combination  of  computational  fluid  dynamics  with  numerical 
optimization  methods.  Some  of  the  earliest  studies  of  such  an  approach  were  made  by  Hicks  and  Henne  [1,2]. 
The  principal  obstacle  was  the  large  computational  cost  of  determining  the  sensitivity  of  the  cost  function 
to  variations  of  the  design  parameters  by  repeated  calculation  of  the  flow.  Another  way  to  approach  the 
problem  is  to  formulate  aerodynamic  shape  design  within  the  framework  of  the  mathematical  theory  for  the 
control  of  systems  governed  by  partial  differential  equations  [3].  In  this  view  the  wing  is  regarded  as  a  device 
to  produce  lift  by  controlling  the  flow,  and  its  design  is  regarded  as  a  problem  in  the  optimal  control  of  the 
flow  equations  by  changing  the  shape  of  the  boundary.  If  the  boundary  shape  is  regarded  as  arbitrary  within 
some  requirements  of  smoothness,  then  the  full  generality  of  shapes  cannot  be  defined  with  a  finite  number 
of  parameters,  and  one  must  use  the  concept  of  the  Frechet  derivative  of  the  cost  with  respect  to  a  function. 
Clearly  such  a  derivative  cannot  be  determined  directly  by  separate  variation  of  each  design  parameter, 
because  there  are  now  an  infinite  number  of  these. 
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Using  techniques  of  control  theory,  however,  the  gradient  can  be  determined  indirectly  by  solving  an 
adjoint  equation  which  has  coefficients  determined  by  the  solution  of  the  flow  equations.  This  directly  cor¬ 
responds  to  the  gradient  technique  for  trajectory  optimization  pioneered  by  Bryson  [4].  The  cost  of  solving 
the  adjoint  equation  is  comparable  to  the  cost  of  solving  the  flow  equations,  with  the  consequence  that  the 
gradient  with  respect  to  an  arbitrarily  large  number  of  parameters  can  be  calculated  with  roughly  the  same 
computational  cost  as  two  flow  solutions.  Once  the  gradient  has  been  calculated,  a  descent  method  can  be 
used  to  determine  a  shape  change  which  will  make  an  improvement  in  the  design.  The  gradient  can  then 
be  recalculated,  and  the  whole  process  can  be  repeated  until  the  design  converges  to  an  optimum  solution, 
usually  within  50  to  100  cycles.  The  fast  calculation  of  the  gradients  makes  optimization  computationally  fea¬ 
sible  even  for  designs  in  three-dimensional  viscous  flow.  There  is  a  possibility  that  the  descent  method  could 
converge  to  a  local  minimum  rather  than  the  global  optimum  solution.  In  practice  this  has  not  proved  a  diffi¬ 
culty,  provided  care  is  taken  in  the  choice  of  a  cost  function  which  properly  reflects  the  design  requirements. 
Conceptually,  with  this  approach  the  problem  is  viewed  as  infinitely  dimensional,  with  the  control  being  the 
shape  of  the  bounding  surface.  Eventually  the  equations  must  be  discretized  for  a  numerical  implementation 
of  the  method.  For  this  purpose  the  flow  and  adjoint  equations  may  either  be  separately  discretized  from 
their  representations  as  differential  equations,  or,  alternatively,  the  flow  equations  may  be  discretized  first, 
and  the  discrete  adjoint  equations  then  derived  directly  from  the  discrete  flow  equations. 

The  effectiveness  of  optimization  as  a  tool  for  aerodynamic  design  also  depends  crucially  on  the  proper 
choice  of  cost  functions  and  constraints.  One  popular  approach  is  to  define  a  target  pressure  distribution, 
and  then  solve  the  inverse  problem  of  finding  the  shape  that  will  produce  that  pressure  distribution.  Since 
such  a  shape  does  not  necessarily  exist,  direct  inverse  methods  may  be  ill-posed.  The  problem  of  designing  a 
two-dimensional  profile  to  attain  a  desired  pressure  distribution  was  studied  by  Lighthill,  who  solved  it  for 
the  case  of  incompressible  flow  with  a  conformal  mapping  of  the  profile  to  a  unit  circle  (5].  The  speed  over 


the  profile  is 


q  =  ^<t>\, 


where  <f>  is  the  potential  which  is  known  for  incompressible  flow  and  h  is  the  modulus  of  the  mapping  function. 
The  surface  value  of  h  can  be  obtained  by  setting  q  =  qd,  where  qd  is  the  desired  speed,  and  since  the  mapping 
function  is  analytic,  it  is  uniquely  determined  by  the  value  of  h  on  the  boundary.  A  solution  exists  for  a 
given  speed  <7oo  at  infinity  only  if 

^jqd.0  =  qoo, 


and  there  are  additional  constraints  on  q  if  the  profile  is  required  to  be  closed. 

The  difficulty  that  the  target  pressure  may  be  unattainable  may  be  circumvented  by  treating  the  inverse 
problem  as  a  special  case  of  the  optimization  problem,  with  a  cost  function  which  measures  the  error  in  the 
solution  of  the  inverse  problem.  For  example,  if  pd  is  the  desired  surface  pressure,  one  may  take  the  cost 
function  to  be  an  integral  over  the  the  body  surface  of  the  square  of  the  pressure  error, 


1  =  ^  Jjp  ~  Pd)2 dB, 

or  possibly  a  more  general  Sobolev  norm  of  the  pressure  error.  This  has  the  advantage  of  converting  a 
possibly  ill  posed  problem  into  a  well  posed  one.  It  has  the  disadvantage  that  it  incurs  the  computational 
costs  associated  with  optimization  procedures. 

The  inverse  problem  still  leaves  the  definition  of  an  appropriate  pressure  architecture  to  the  designer. 
One  may  prefer  to  directly  improve  suitable  performance  parameters,  for  example,  to  minimize  the  drag  at  a 
given  lift  and  Mach  number.  In  this  case  it  is  important  to  introduce  appropriate  constraints.  For  example,  \ 
if  the  span  is  not  fixed  the  vortex  drag  can  be  made  arbitrarily  small  by  sufficiently  increasing  the  span.  In 
practice,  a  useful  approach  is  to  fix  the  planform,  and  optimize  the  wing  sections  subject  to  constraints  on 
minimum  thickness. 


2  Formulation  of  the  Design  Problem  as  a  Control  Problem 

The  simplest  approach  to  optimization  is  to  define  the  geometry  through  a  set  of  design  parameters,  which 
may,  for  example,  be  the  weights  o,  applied  to  a  set  of  shape  functions  bi{x)  so  that  the  shape  is  represented 


CIO  »  1 1  ^ 

fix)  =  5 Zaibi(x )• 

Then  a  cost  function  I  is  selected  which  might,  for  example,  be  the  drag  coefficient  or  the  lift  to  drag  ratio, 
and  I  is  regarded  as  a  function  of  the  parameters  at-  The  sensitivities  may  now  be  estimated  by  making 
a  small  variation  Sa,  in  each  design  parameter  in  turn  and  recalculating  the  flow  to  obtain  the  change  in  I. 

Then  dl_  _  Ijoj  +  Sotj)  -  Ija i) 

doci  Scti 

The  gradient  vector  may  now  be  used  to  determine  a  direction  of  improvement.  The  simplest  procedure 
is  to  make  a  step  in  the  negative  gradient  direction  by  setting 

an+i  =  an  —  A  Ja, 


so  that  to  first  order  T  ^jT  ^ 

I +  81  =  1-  -L-8a  =  1-  . 

da  da  da 

More  sophisticated  search  procedures  may  be  used  such  as  quasi-Newton  methods,  which  attempt  to  estimate 
the  second  derivative  of  the  cost  function  from  changes  in  the  gradient  in  successive  optimization 

steps.  These  methods  also  generally  introduce  line  searches  to  find  the  minimum  in  the  search  direction  which 
is  defined  at  each  step.  The  main  disadvantage  of  this  approach  is  the  need  for  a  number  of  flow  calculations 
proportional  to  the  number  of  design  variables  to  estimate  the  gradient.  The  computational  costs  can  thus 
become  prohibitive  as  the  number  of  design  variables  is  increased. 

Using  techniques  of  control  theory,  however,  the  gradient  can  be  determined  indirectly  by  solving  an 
adjoint  equation  which  has  coefficients  defined  by  the  solution  of  the  flow  equations.  The  cost  of  solving  the 
adjoint  equation  is  comparable  to  that  of  solving  the  flow  equations.  Thus  the  gradient  can  be  determined 
with  roughly  the  computational  costs  of  two  flow  solutions,  independently  of  the  number  of  design  variables, 
which  may  be  infinite  if  the  boundary  is  regarded  as  a  free  surface.  The  underlying  concepts  are  clarified  by 
the  following  abstract  description  of  the  adjoint  method. 

For  flow  about  an  airfoil  or  wing,  the  aerodynamic  properties  which  define  the  cost  function  are  functions 
of  the  flow-field  variables  (w)  and  the  physical  location  of  the  boundary,  which  may  be  represented  by  the 
function  T,  say.  Then 

I  =  I  {w,F) , 

and  a  change  in  T  results  in  a  change 


in  the  cost  function.  Here,  the  subscripts  1  and  JJ  are  used  to  distinguish  the  contributions  due  to  the 
variation  8w  in  the  flow  solution  from  the  change  associated  directly  with  the  modification  bT  in  the  shape. 
This  notation  assists  in  grouping  the  numerous  terms  that  arise  during  the  derivation  of  the  full  Navier 
Stokes  adjoint  operator,  outlined  later,  so  that  the  basic  structure  of  the  approach  as  it  is  sketched  in  the 
present  section  can  easily  be  recognized. 

Suppose  that  the  governing  equation  R  which  expresses  the  dependence  of  w  and  T  within  the  flowfield 
domain  D  can  be  written  as 

R  (tu,  T)  =  0.  (2) 


Then  Sw  is  determined  from  the  equation 


to+[§pj  ^  =  0- 

dw  j  J  II 


(3) 


Since  the  variation  SR  is  zero,  it  can  be  multiplied  by  a  Lagrange  Multiplier  ip  and  subtracted  from  the 
variation  SI  without  changing  the  result.  Thus  equation  (1)  can  be  replaced  by 


S  dw  W 


ST. 


(4) 


Choosing  rp  to  satisfy  the  adjoint  equation 


—]T it  -  — 

dw  dw 


(5) 


the  first  term  is  eliminated,  and  we  find  that 


SI  =  GST, 


(6) 


where 


The  advantage  is  that  (6)  is  independent  of  Sw ,  with  the  result  that  the  gradient  of  I  with  respect  to  an 
arbitrary  number  of  design  variables  can  be  determined  without  the  need  for  additional  flow-field  evaluations. 
In  the  case  that  (2)  is  a  partial  differential  equation,  the  adjoint  equation  (5)  is  also  a  partial  differential 
equation  and  determination  of  the  appropriate  boundary  conditions  requires  careful  mathematical  treatment. 

In  reference  [6]  Jameson  derived  the  adjoint  equations  for  transonic  flows  modeled  by  both  the  potential 
flow  equation  and  the  Euler  equations.  The  theory  was  developed  in  terms  of  partial  differential  equations, 
leading  to  an  adjoint  partial  differential  equation.  In  order  to  obtain  numerical  solutions  both  the  flow  and 
the  adjoint  equations  must  be  discretized.  The  control  theory  might  be  applied  directly  to  the  discrete  flow 
equations  which  result  from  the  numerical  approximation  of  the  flow  equations  by  finite  element,  finite 
volume  or  finite  difference  procedures.  This  leads  directly  to  a  set  of  discrete  adjoint  equations  with  a  matrix 
which  is  the  transpose  of  the  Jacobian  matrix  of  the  full  set  of  discrete  nonlinear  flow  equations.  On  a 
three-dimensional  mesh  with  indices  i7j,k  the  individual  adjoint  equations  may  be  derived  by  collecting 
together  all  the  terms  multiplied  by  the  variation  Switjyk  of  the  discrete  flow  variable  The  resulting 

discrete  adjoint  equations  represent  a  possible  discretization  of  the  adjoint  partial  differential  equation.  If 
these  equations  are  solved  exactly  they  can  provide  an  exact  gradient  of  the  inexact  cost  function  which 
results  from  the  discretization  of  the  flow  equations.  The  discrete  adjoint  equations  derived  directly  from  the 
discrete  flow  equations  become  very  complicated  when  the  flow  equations  are  discretized  with  higher  order 
upwind  biased  schemes  using  flux  limiters.  On  the  other  hand  any  consistent  discretization  of  the  adjoint 
partial  differential  equation  will  yield  the  exact  gradient  in  the  limit  as  the  mesh  is  refined.  The  trade-off 
between  the  complexity  of  the  adjoint  discretization,  the  accuracy  of  the  resulting  estimate  of  the  gradient, 
and  its  impact  on  the  computational  cost  to  approach  an  optimum  solution  is  a  subject  of  agoing  research. 

The  true  optimum  shape  belongs  to  an  infinitely  dimensional  space  of  design  parameter^  One  motivation 
for  developing  the  theory  for  the  partial  differential  equations  of  the  flow  is  to  provide  an  indication  in 
principle  of  how  such  a  solution  could  be  approached  if  sufficient  computational  resources  were  available. 
Another  motivation  is  that  it  highlights  the  possibility  of  generating  ill  posed  formulations  of  the  problem. 
For  example,  if  one  attempts  to  calculate  the  sensitivity  of  the  pressure  at  a  particular  location  to  changes 
in  the  boundary  shape,  there  is  the  possibility  that  a  shape  modification  could  cause  a  shock  wave  to  pass 
over  that  location.  Then  the  sensitivity  could  become  unbounded.  The  movement  of  the  shock,  however, 
is  continuous  as  the  shape  changes.  Therefore  a  quantity  such  as  the  drag  coefficient,  which  is  determined 
by  integrating  the  pressure  over  the  surface,  also  depends  continuously  on  the  shape.  The  adjoint  equation 
allows  the  sensitivity  of  the  drag  coefficient  to  be  determined  without  the  explicit  evaluation  of  pressure 
sensitivities  which  would  be  ill  posed. 

The  discrete  adjoint  equations,  whether  they  are  derived  directly  or  by  discretization  of  the  adjoint 
partial  differential  equation,  are  linear.  Therefore  they  could  be  solved  by  direct  numerical  inversion.  In 
three-dimensional  problems  on  a  mesh  with,  say,  n  intervals  in  each  coordinate  direction,  the  number  of 
unknowns  is  proportional  to  n3  and  the  bandwidth  to  n2.  The  complexity  of  direct  inversion  is  proportional 
to  the  number  of  unknowns  multiplied  by  the  square  of  the  bandwidth,  resulting  in  a  complexity  proportional 
to  n7.  The  cost  of  direct  inversion  can  thus  become  prohibitive  as  the  mesh  is  refined,  and  it  becomes  more 
efficient  to  use  iterative  solution  methods.  Moreover,  because  of  the  similarity  of  the  adjoint  equations  to 
the  flow  equations,  the  same  iterative  methods  which  have  been  proved  to  be  efficient  for  the  solution  of  the 
flow  equations  are  efficient  for  the  solution  of  the  adjoint  equations. 


Studies  of  the  use  of  control  theory  for  optimum  shape  design  of  systems  governed  by  elliptic  equations 
were  initiated  by  Pironneau  [7].  The  control  theory  approach  to  optimal  aerodynamic  design  was  first  applied 
to  transonic  flow  by  Jameson  [8,6,9-12],  He  formulated  the  method  for  inviscid  compressible  flows  with 
shock  waves  governed  by  both  the  potential  flow  and  the  Euler  equations  [6],  Numerical  results  showing  the 
method  to  be  extremely  effective  for  the  design  of  airfoils  in  transonic  potential  flow  were  presented  in  [13], 
and  for  three-dimensional  wing  design  using  the  Euler  equations  in  [14],  More  recently  the  method  has  been 
employed  for  the  shape  design  of  complex  aircraft  configurations  [15,16],  using  a  grid  perturbation  approach 
to  accommodate  the  geometry  modifications.  The  method  has  been  used  to  support  the  aerodynamic  design 
studies  of  several  industrial  projects,  including  the  Beech  Premier  and  the  McDonnell  Douglas  MDXX  and 
Blended  Wing-Body  projects.  The  application  to  the  MDXX  is  described  in  [10].  The  experience  gained  in 
these  industrial  applications  made  it  clear  that  the  viscous  effects  cannot  be  ignored  in  transonic  wing  design, 
and  the  method  has  therefore  been  extended  to  treat  the  Reynolds  Averaged  Navier-Stokes  equations  [12]. 
Adjoint  methods  have  also  been  the  subject  of  studies  by  a  number  of  other  authors,  including  Baysal  and 
Eleshaky  [17],  Huan  and  Modi  [18],  Desai  and  Ito  [19],  Anderson  and  Venkatakrishnan  [20],  and  Peraire 
and  Elliot  [21].  Ta’asan,  Kuruvila  and  Salas  [22],  who  have  implemented  a  one  shot  approach  in  which 
the  constraint  represented  by  the  flow  equations  is  only  required  to  be  satisfied  by  the  final  converged 
solution.  In  their  work,  computational  costs  are  also  reduced  by  applying  multigrid  techniques  to  the  geometry 
modifications  as  well  as  the  solution  of  the  flow  and  adjoint  equations. 

The  next  sections  discuss  the  application  of  the  method  to  automatic  wing  design  with  the  flow  modeled 
by  the  three-dimensional  Euler  and  Navier-Stokes  equations.  The  computational  costs  are  low  enough  that 
it  has  proved  possible  to  determine  optimum  wing  designs  using  the  Euler  equations  in  a  few  hours  on 
workstations  such  as  the  IBM590  or  the  Silicon  Graphics  Octane. 


3  The  Navier-Stokes  Equations 


For  the  derivations  that  follow,  it  is  convenient  to  use  Cartesian  coordinates  (xi,x2,x3)  and  to  adopt  the 
convention  of  indicial  notation  where  a  repeated  index  V  implies  summation  over  i  =  1  to  3.  The  three- 
dimensional  Navier-Stokes  equations  then  take  the  form 


dw  df± 
dt  dxi 


dfyj 

dXi 


in  V, 


(7)  . 


where  the  state  vector  w,  inviscid  flux  vector  /  and  viscous  flux  vector  /„  are  described  respectively  by 
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r  \ 

pUi 
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pui 

pUiUi  -1-  p5ii 

CTijSji 
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pUiU2  +  pSi2 
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&  ij  2  ► 

PU  3 

pUiUz  +  pSi3 

®ij  3 

,  ujaij  +  J 

[pE  ) 

prnH 

In  these  definitions,  p  is  the  density,  u\ ,  uo,  tt3  are  the  Cartesian  velocity  components,  E  is  the  total  energy 
and  6ij  is  the  Kronecker  delta  function.  The  pressure  is  determined  by  the  equation  of  state 


P  = 


and  the  stagnation  enthalpy  is  given  by 


H  =  E+~, 
P 


(9) 


where  7  is  the  ratio  of  the  specific  heats.  The  viscous  stresses  may  be  written  as 

( dm  ,dUj\,  vr  duk 

c«  “  "  +  tel)  + 

where  p  and  A  are  the  first  and  second  coefficients  of  viscosity.  The  coefficient  of  thermal  conductivity  and 
the  temperature  are  computed  as 

l-EeHt-JL 

Pr1  Rp 


(10) 


where  Pr  is  the  Prandtl  number,  cp  is  the  specific  heat  at  constant  pressure,  and  R  is  the  gas  constant. 

For  discussion  of  real  applications  using  a  discretization  on  a  body  conforming  structured  mesh,  it  is  also 
useful  to  consider  a  transformation  to  the  computational  coordinates  (6.6.6)  defined  by  the  metrics 


J  —  det  (K) ,  K^  = 


■%' 

dxj 


The  Navier-Stokes  equations  can  then  be  written  in  computational  space  as 


d(Jw)  d(Fx-  Fm)  =  {n) 

at  d£i 

where  the  inviscid  and  viscous  flux  contributions  are  now  defined  with  respect  to  the  computational  cell 
faces  by  =  S0/j  and  Fvi  =  %fri,  and  the  quantity  %  =  JK~^  represents  the  projection  of  the  6  cell 
face  along  the  xj  axis.  In  obtaining  equation  (11)  we  have  made  use  of  the  property  that 

?pi  =  o  (12) 

which  represents  the  fact  that  the  sum  of  the  face  areas  over  a  closed  volume  is  zero,  as  can  be  readily  verified 
by  a  direct  examination  of  the  metric  terms. 


4  Formulation  of  the  Optimal  Design  Problem  for  the  Navier-Stokes  Equations 


Aerodynamic  optimization  is  based  on  the  determination  of  the  effect  of  shape  modifications  on  some  perfor¬ 
mance  measure  which  depends  on  the  flow.  For  convenience,  the  coordinates  £*  describing  the  fixed  compu¬ 
tational  domain  are  chosen  so  that  each  boundary  conforms  to  a  constant  value  of  one  of  these  coordinates. 
Variations  in  the  shape  then  result  in  corresponding  variations  in  the  mapping  derivatives  defined  by  K{j. 
Suppose  that  the  performance  is  measured  by  a  cost  function 

1  =  J  M{w,S)dBi  +  J  V(w,S)dDt, 

containing  both  boundary  and  field  contributions  where  dB^  and  dD^  are  the  surface  and  volume  elements 
in  the  computational  domain.  In  general,  M  and  V  will  depend  on  both  the  flow  variables  w  and  the  metrics 
S  defining*  the  computational  space.  The  design  problem  is  now  treated  as  a  control  problem  where  the 
boundary  shape  represents  the  control  function,  which  is  chosen  to  minimize  I  subject  to  the  constraints 
defined  by  the  flow  equations  (11).  A  shape  change  produces  a  variation  in  the  flow  solution  6w  and  the 
metrics  5S  which  in  turn  produce  a  variation  in  the  cost  function 

61=  [  6M{w,  S)  dB$  4-  [  6V(w,S)dDz.  (13) 

Jb  Jv 


This  can  be  split  as 


SI  =  61  i  4-  81  a. 


(14) 


with 


8 M  =  [Mw\i  8w  +  SMn , 

6P  =  [Pw]I6w  +  SPII,  (15) 

where  we  continue  to  use  the  subscripts  I  and  H  to  distinguish  between  the  contributions  associated  with 
the  variation  of  the  flow  solution  Sw  and  those  associated  with  the  metric  variations  6S.  Thus  [Mi.]/  and 
[Pw\j  represent  ^7  and  with  the  metrics  fixed,  while  SMn  and  SVu  represent  the  contribution  of  the 
metric  variations  8S  to  8M  and  8V . 


In  the  steady  state,  the  constraint  equation  (53)  specifies  the  variation  of  the  state  vector  Sw  by 


SR=~S(Fi-Fvi)  =  0.  r 


(16) 


Here,  also,  SR,  SFi  and  SFvi  can  be  split  into  contributions  associated  with  Sw  and  SS  using  the  notation 

SR  —  SRj  +  SRjj 
SFi  =  [Fiw},Sw  +  SFm 

SFvi  =  [iViio]/  Sw  +  SFyiij.  (17) 

The  inviscid  contributions  are  easily  evaluated  as 

[*.],  =  Sij^,  6FviH  =  SSijfj. 


The  details  of  the  viscous  contributions  are  complicated  by  the  additional  level  of  derivatives  in  the  stress 
and  heat  flux  terms. 

Multiplying  by  a  co-state  vector  rp,  which  will  play  an  analogous  role  to  the  Lagrange  multiplier  introduced 
in  equation  (4),  and  integrating  over  the  domain  produces 


f  ■tpT^F8{Fi-Fvi)dVi  =  0. 

Jv 


dii 

Assuming  that  x p  is  differentiable  the  terms  with  subscript  I  may  be  integrated  by  parts  to  give 

/,T 


(18) 


/  mVS  (. Fi  -  Fvi)j  dB{  -  f  5  ( Fi  -  Fvi)j  dPc  +  f  yf SRudV €  =  0.  (19) 

JB  JT>  aSi  Jv 


This  equation  results  directly  from  talcing  the  variation  of  the  weak  form  of  the  flow  equations,  where  xp  is 
taken  to  be  an  arbitrary  differentiable  test  function.  Since  the  left  hand  expression  equals  zero,  it  may  be 
subtracted  from  the  variation  in  the  cost  function  (13)  to  give 


SI  =  5I„  -  f  xpTSR„dVi  -  f  [5 Ml  -  n^S  (Fi  -  Fvi)/]  dBi 
Jv  Jb 

+  [sPi  +  ^-6  (Fi  -  F„),]  (20) 

Now,  since  xp  is  an  arbitrary  differenaable  function,  it  may  be  chosen  in  such  a  way  that  SI  no  longer 
depends  explicitly  on  the  variation  of  the  state  vector  Sw.  The  gradient  of  the  cost  function  can  then  be 
evaluated  directly  from  the  metric  variations  without  having  to  recompute  the  variation  Sw  resulting  from 
the  perturbation  of  each  design  variable. 

Comparing  equations  (15)  and  (17),  the  variation  Sw  may  be  eliminated  from  (20)  by  equating  all  field 
terms  with  subscript  UP  to  produce  a  differential  adjoint  system  governing  xp 


r\i  r 

-qT"  [Fiw  -  Fmu,]/  +  [pw\i  =  0  in  V. 


(21) 


The  corresponding  adjoint  boundary  condition  is  produced  by  equating  the  subscript  “J”  boundary  terms 
in  equation  (20)  to  produce 

rc«V,T  [-Fi™  -  Fviw]j  =  \MJ\i  on  B.  (22) 

The  remaining  terms  from  equation  (20)  then  yield  a  simplified  expression  for  the  variation  of  the  cost 
function  which  defines  the  gradient 

SI  =  Sin  +  /  VSRndV^,  (23) 

Jv 

which  consists  purely  of  the  terms  containing  variations  in  the  metrics  with  the  flow  solution  fixed.  Hence 
an  explicit  formula  for  the  gradient  can  be  derived  once  the  relationship  between  mesh  perturbations  and 
shape  variations  is  defined. 


Comparing  equations  (15)  and  (17),  the  variation  5w  may  be  eliminated  from  (20)  by  equating  all  field 
terms  with  subscript  to  produce  a  differential  adjoint  system  governing  ip 

[Fiw  -  Fviw]j  +VW  =  0  mV.  (24) 

The  corresponding  adjoint  boundary  condition  is  produced  by  equating  the  subscript  boundary  terms 
in  equation  (20)  to  produce 

mipT  [Fil0  -  Fviw}j  =  MW  on  B.  (25) 

The  remaining  terms  from  equation  (20)  then  yield  a  simplified  expression  for  the  variation  of  the  cost 
function  which  defines  the  gradient 

SI  =  J  {SMn  -  riiipT  [<SFi  -  <5F„i]  a}  dB 4  -\J  j SVn  +  [6Fj  -  5F„j]  //|  dV$.  (26) 

The  details  of  the  formula  for  the  gradient  depend  on  the  way  in  which  the  boundary  shape  is  parameterized  as 
a  function  of  the  design  variables,  and  the  way  in  which  the  mesh  is  deformed  as  the  boundary  is  modified. 
Using  the  relationship  between  the  mesh  deformation  and  the  surface  modification,  the  field  integral  is 
reduced  to  a  surface  integral  by  integrating  along  the  coordinate  lines  emanating  from  the  surface.  Thus  the 
expression  for  SI  is  finally  reduced  to  the  form  of  equation  (6) 

61=  [  QSTdBz 

Jb 

where  T  represents  the  design  variables,  and  Q  is  the  gradient,  which  is  a  function  defined  over  the  boundary 
surface. 

The  boundary  conditions  satisfied  by  the  flow  equations  restrict  the  form  of  the  left  hand  side  of  the 
adjoint  boundary  condition  (25).  Consequently,  the  boundary  contribution  to  the  cost  function  M  carinot 
be  specified  arbitrarily.  Instead,  it  must  be  chosen  from  the  class  of  functions  which  allow  cancellation  of 
all  terms  containing  Sw  in  the  boundary  integral  of  equation  (20).  On  the  other  hand,  there  is  no  such 
restriction  on  the  specification  of  the  field  contribution  to  the  cost  function  V ,  since  these  terms  may  always 
be  absorbed  into  the  adjoint  field  equation  (24)  as  source  terms. 

It  is  convenient  to  develop  the  inviscid  and  viscous  contributions  to  the  adjoint  equations  separately.  Also, 
for  simplicity,  it  will  be  assumed  that  the  portion  of  the  boundary  that  undergoes  shape  modifications  is 
restricted  to  the  coordinate  surface  £2  =  0.  Then  equations  (20)  a  id  (22)  may  be  simplified  by  incorporating 
the  conditions 

Tt\  713  —  0,  7t»2  1,  dB £  d^\d^ 3? 

so  that  only  the  variations  SF2  and  6Fv2  need  to  be  considered  at  the  wall  boundary. 


5  Derivation  of  the  Inviscid  Adjoint  Terms 


The  inviscid  contributions  have  been  previously  derived  in  [13,23]  but  are  included  here  for  completeness. 
Taking  the  transpose  of  equation  (21),  the  inviscid  adjoint  equation  may  be  written  as 

Cff^=0  in  V,  (27) 

where  the  inviscid  Jacobian  matrices  in  the  transformed  space  are  given  by 


C* 


Sij 


The  transformed  velocity  components  have  the  fomv 


and  the  condition  that  there  is  no  flow  through  the  wall  boundary  at  &  =  0  is  equivalent  to 


U2=  0, 


so  that 

SU2=0 

when  the  boundary  shape  is  modified.  Consequently  the  variation  of  the  inviscid  flux  at  the  boundary  reduces 
to 


SF2  =  Sp  < 
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(28) 


Since  SF2  depends  only  on  the  pressure,  it  is  now  clear  that  the  performance  measure  on  the  boundary 
_A/f  (u;,  S)  may  only  be  a  function  of  the  pressure  and  metric  terms.  Otherwise,  complete  cancellation  of  the 
terms  containing  5w  in  the  boundary  integral  would  be  impossible.  One  may,  for  example,  include  arbitrary 
measures  of  the  forces  and  moments  in  the  cost  function,  since  these  are  functions  of  the  surface  pressure. 
In  order  to  design  a  shape  which  will  lead  to  a  desired  pressure  distribution,  a  natural  choice  is  to  set 


dS 


where  pd  is  the  desired  surface  pressure,  and  the  integral  is  evaluated  over  the  actual  surface  area.  In  the 
computational  domain  this  is  transformed  to 


=  \JJB  (p-Pd)2|S2|<^3, 


where  the  quantity  _ 

13,1  =  y/S^STj 

denotes  the  face  area  corresponding  to  a  unit  element  of  face  area  in  the  computational  domain.  Now,  to 
cancel  the  dependence  of  the  boundary  integral  on  Sp,  the  adjoint  boundary  condition  reduces  to 

tfjrij  =  p  —  Pd  (29) 


where  n3  are  the  components  of  the  surface  normal 


rii  = 


_  02  j 


i&r 


This  amounts  to  a  transpiration  boundary  condition  on  the  co-state  variables  corresponding  to  the  momen¬ 
tum  components.  Note  that  it  imposes  no  restriction  on  the  tangential  component  of  *0  &t  the  boundary. 

In  the  presence  of  shock  waves,  neither  p  nor  pd  are  necessarily  continuous  at  the  surface.  The  boundary 
condition  is  then  in  conflict  with  the  assumption  that  ip  is  differentiable.  This  difficulty  can  be  circumvented 
by  the  use  of  a  smoothed  boundary  condition  [23]. 

6  Derivation  of  the  Viscous  Adjoint  Terms 

In  computational  coordinates,  the  viscous  terms  in  the  N avier— Stokes  equations  have  the  form 


Computing  the  variation  Sw  resulting  from  a  shape  modification  of  the  boundary,  introducing  a  co-state 
vector  ip  and  integrating  by  parts  following  the  steps  outlined  by  equations  (16)  to  (19)  produces 

J  ipT  [SSijfvj  +  S2jSfvj)  dB^  —  (5SijfVj  +  SijSfvj)  dD^, 

where  the  shape  modification  is  restricted  to  the  coordinate  surface  (2  =  0  so  that  n\  =  713  =  0,  and  ri2  =  1. 
Furthermore,  it  is  assumed  that  the  boundary  contributions  at  the  far  field  may  either  be  neglected  or  else 
eliminated  by  a  proper  choice  of  boundary  conditions  as  previously  shown  for  the  inviscid  case  [13,23]. 

The  viscous  terms  will  be  derived  under  the  assumption  that  the  viscosity  and  heat  conduction  coefficients 
/x  and  k  are  essentially  independent  of  the  flow,  and  that  their  variations  may  be  neglected.  This  simplification 
has  been  successfully  used  for  may  aerodynamic  problems  of  interest.  In  the  case  of  some  turbulent  flows, 
there  is  the  possibility  that  the  flow  variations  could  result  in  significant  changes  in  the  turbulent  viscosity, 
and  it  may  then  be  necessary  to  account  for  its  variation  in  the  calculation. 


Transformation  to  Primitive  Variables 

The  derivation  of  the  viscous  adjoint  terms  is  simplified  by  transforming  to  the  primitive  variables 

wT  =  (p,Ui,u2jU3,p)t1 

because  the  viscous  stresses  depend  on  the  velocity  derivatives  ,  while  the  heat  flux  can  be  expressed  as 


K 


dxi  \p 


where  k  =  ^  =  p^zrpj-  The  relationship  between  the  conservative  and  primitive  variations  is  defined  by 
the  expressions 

Sw  =  MSw ,  Sw  =  M~1Sw 

which  make  use  of  the  transformation  matrices  M  =  §f  and  M^1  =  ff .  These  matrices  are  provided  in 
transposed  form  for  future  convenience 
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The  conservative  and  primitive  adjoint  operators  L  and  L  corresponding  to  the  variations  6w  and  Sw  are 
then  related  by 


f  5wTLip  dV^  =  f  SwT  Lip  dV^ 
Jv  Jv 


with 

L  =  MTLi 


so  that  after  determining  the  primitive  adjoint  operator  by  direct  evaluation  of  the  viscous  portion  of  (21), 
the  conservative  operator  may  be  obtained  by  the  transformation  L  =  L.  Since  the  continuity  equation 

contains  no  viscous  terms,  it  makes  no  contribution  to  the  viscous  adjoint  system.  Therefore,  the  derivation 
proceeds  by  first  examining  the  adjoint  operators  arising  from  the  momentum  equations. 


Contributions  from  the  Momentum  Equations 

In  order  to  make  use  of  the  summation  convention,  it  is  convenient  to  set  ijjj+ i  =  <f>j  for  j  =  1,2,3.  Then 
the  contribution  from  the  momentum  equations  is 

J  4>k  (SS2jakj  +  S2j8akj)  (SSym  +  Si:jSakj)  dV^.  (30) 

The  velocity  derivatives  in  the  viscous  stresses  can  be  expressed  as 

dui  __  dui  d£i  _  Sjj  du^ 

dxj  ~  d£i  dxj  J 

with  corresponding  variations 


dui  rSyl  d  r^Uil  r  (Sij\ 


The  variations  in  the  stresses  are  then 

Sakj  =  (m  [^j- ^Suk  +  m;Sui\  +  A  lSjk f5r4i'<5Um]  } / 

+  {m  [s  (%*)  + 6  (^fc)  §&  ] +  A  I5’*5  )  %r] }//  ■ 

As  before  only  those  terms  with  subscript  I,  which  contain  variations  of  the  flow  variables,  need  be  considered 
further  in  deriving  the  adjoint  operator.  The  field  contributions  that  contain  8Ui  in  equation  (30)  appear  as 

-  L  {*  (H* +  H*) 

This  may  be  integrated  by  parts  to  yield 

+  IvSujWi{SlkSiijW)dD* 

where  the  boundary  integral  has  been  eliminated  by  noting  that  Sui  =  0  on  the  solid  boundary.  By  exchanging 
indices,  the  field  integrals  may  be  combined  to  produce 


Sij  d<t>k  Sjk  d<j)j  \ 

“86  +  T  86  > 


im  S$m  \ 

TWI™1' 


which  is  further  simplified  by  transforming  the  inner  derivatives  back  to  Cartesian  coordinates 

The  boundary  contributions  that  contain  8u{  in  equation  (30)  may  be  simplified  using  the  fact  that 


-£r5ui=0  if  1  =  1,3 


on  the  boundary  B  so  that  they  become 


Together,  (31)  and  (32)  comprise  the  field  and  boundary  contributions  of  the  momentum  equations  to  the 
viscous  adjoint  operator  in  primitive  variables. 


Contributions  from  the  Energy  Equation 

In  order  to  derive  the  contribution  of  the  energy  equation  to  the  viscous  adjoint  terms  it  is  convenient  to  set 

+*  =  0>  Q*  =  Uiai*  +  (p)  ’ 

where  the  temperature  has  been  written  in  terms  of  pressure  and  density  using  (10).  The  contribution  from 
the  energy  equation  can  then  be  written  as 

J  e  ( sSijQj  +  s2jSQj)  dBt~  (5SijQj +  svsQi)  dVt-  (33) 

The  field  contributions  that  contain  8uu8p ,  and  8p  in  equation  (33)  appear  as 

■  L  w<s‘iSQitDi  -  -  Lw,s'j  i6*™ + ukSni  (%  ~  ??)) <n>t' 

The  term  involving  Scrkj  may  ^>e  integrated  by  parts  to  produce 

LsukwiSij  Mufc^+u,S) 

where  the  conditions  it<  =  8m  =  0  are  used  to  eliminate  the  boundary  integral  on  B.  Notice  that  the  other 
term  in  (34)  that  involves  8uk  need  not  be  integrated  by  parts  and  is  merely  carried  on  as 


(34) 


(35) 


r  ae 


(36) 


The  terms  in  expression  (34)  that  involve  8p  and  Sp  may  also  be  integrated  by  parts  to  produce  both  a 
field  and  a  boundary  integral.  The  field  integral  becomes 

which  may  be  simplified  by  transform?  ug^he  inner  derivative  to  Cartesian  coordinates 

I  (6-E-Ps-£)  ±  (SljK ~)d v,. 

Jt>\P.  P  P  )  \  / 

The  boundary  integral  becomes 

I JSJL  _  £<2\ 

Jb  \  p  p  p  J  3  dZ' 

This  can  be  simplified  by  transforming  the  inner  derivative  to  Cartesian  coordinates 

fJk.Ek) 

Jb  \P  P  P  )  3  dxi 


(37) 


(38) 


(39) 


and  identifying  the  normal  derivative  at  the  wall 


±  =  c  JL 

dn  2j  dxj  ’ 


OT.i  (*-2&), 


(40) 


and  the  variation  in  temperature 


to  produce  the  boundary  contribution 


(41) 


/.<*<•  F 

This  term  vanishes  if  T  is  constant  on  the  wall  but  persists  if  the  wall  is  adiabatic. 

There  is  also  a  boundary  contribution  left  over  from  the  first  integration  by  parts  (33)  which  has  the 


form 


(42) 


where 


Qi 


=  k 


dT 

dXj' 


since  u,  =  0.  Notice  that  for  future  convenience  in  discussing  the  adjoint  boundary  conditions  resulting  from 
the  energy  equation,  both  the  Sw  and  SS  terms  corresponding  to  subscript  classes  I  and  U  are  considered 
simultaneously.  If  the  wall  is  adiabatic 


so  that  using  (40), 

S  ( S2jQj )  =  o, 

and  both  the  6w  and  6S  boundary  contributions  vanish. 

On  the  other  hand,  if  T  is  constant  =  0  for  l  =  1,3,  so  that 


Qi 


Thus,  the  boundary  integral  (42)  becomes 


Vi 

J  d& 


6T  +  6 


(43) 


Therefore,  for  constant  T,  the  first  term  corresponding  to  variations  in  the  flow  field  contributes  to  the 
adjoint  boundary  operator  and  the  second  set  of  terms  corresponding  to  metric  variations  contribute  to  the 
cost  function  gradient. 

All  together,  the  contributions  from  the  energy  equation  to  the  viscous  adjoint  operator  are  the  three 
field  terms  (35),  (36)  and  (37),  and  either  of  two  boundary  contributions  ( 41)  or  (  43),  depending  on  whether 
the  wall  is  adiabatic  or  has  constant  temperature. 


The  Viscous  Adjoint  Field  Operator 

Collecting  together  the  contributions  from  the  momentum  and  energy  equations,  the  viscous  adjoint  operator 
in  primitive  variables  can  be  expressed  as 


for  i  =  1,2,3 


The  conservative  viscous  adjoint  operator  may  now  be  obtained  by  the  transformation 


L  -  M~lTL. 


7  Viscous  Adjoint  Boundary  Conditions 

It  was  recognized  in  Section  4  that  the  boundary  conditions  satisfied  by  the  flow  equations  restnct  the  form 
of  the  performance  measure  that  may  be  chosen  for  the  cost  function.  There  must  be  a  direct  correspondence 
between  the  flow  variables  for  which  variations  appear  in  the  variation  of  the  cost  function,  and  those  variables 
for  which  variations  appear  in  the  boundary  terms  arising  during  the  derivation  of  the  adjoint  field  equations. 
Otherwise  it  would  be  impossible  to  eliminate  the  dependence  of  SI  on  Sw  through  proper  specification  of 
the  adjoint  boundary  condition.  As  in  the  derivation  of  the  field  equations,  it  proves  convenient  to  consider 
the  contributions  from  the  momentum  equations  and  the  energy  equation  separately. 


Boundary  Conditions  Arising  from  the  Momentum  Equations 

The  boundary  term  that  arises  from  the  momentum  equations  including  both  the  Sw  and  SS  components 
(30)  takes  the  form 

J  <t>k&  {SljVkj)  dB^. 

Replacing  the  metric  term  with  the  corresponding  local  face  area  S2  and  unit  normal  rij  defined  by 

\Sa\  —  y/S^STj,  nj  = 


then  leads  to 


/  4>k&  (l^l  nj&kj)  dB$. 

Jb 


Defining  the  components  of  the  surface  stress  as 


Tfc  —  TijCFkj 

and  the  physical  surface  element 

dS  =  |S2| 

the  integral  may  then  be  split  into  two  components 


f  <t>kTk  |<5jS,2|  dB$  +  [  <f>k$TkdS, 
Jb  Jb 


where  onlv  the  second  term  contains  variations  in  the  flow  variables  and  must  consequently  cancel  the  Sw 
terms  arising  in  the  cost  function.  The  first  term  will  appear  in  the  expression  for  the  gradient. 

A  general  expression  for  the  cost  function  that  allows  cancellation  with  terms  containing  St k  has  the  form 


corresponding  to  a  variation 


r  =  f  A f(r)dS, 

Jb 

~LZSTtds’ 


for  which  cancellation  is  achieved  by  the  adjoint  boundary  condition 

dN 

Natural  choices  for  M  arise  from  force  optimization  and  as  measures  of  the  deviation  of  the  surface  stresses 
from  desired  target  values. 

For  viscous  force  optimization,  the  cost  function  should  measure  friction  drag.  The  friction  force  in  the 
direction  is 


(T%jdSj  — 


S2j  & ij  dB^ 


so  that  the  force  in  a  direction  with  cosines  qt  has  the  form 

Cqf  =  /  <7i S2j<7ijdBc. 

JB 

Expressed  in  terms  of  the  surface  stress  this  corresponds  to 

Cqf  =  /  q^dS , 

JB 

so  that  basing  the  cost  function  (45)  on  this  quantity  gives 

M  =  qin. 

Cancellation  with  the  flow  variation  terms  in  equation  (44)  therefore  mandates  the  adjoint  boundary  condi¬ 
tion 

<t>k  =  ft*. 

Note  that  this  choice  of  boundary  condition  also  eliminates  the  first  term  in  equation  (44)  so  that  it  need 
not  be  included  in  the  gradient  calculation. 

In  the  inverse  design  case,  where  the  cost  function  is  intended  to  measure  the  deviation  of  the  surface 
stresses  from  some  desired  target  values,  a  suitable  definition  is 

M{t)  =  ^aik  in  -  T di)  (jk  -  Tdk) . 

where  rd  is  the  desired  surface  stress,  including  the  contribution  of  the  pressure,  and  the  coefficients  alk 
define  a  weighting  matrix.  For  cancellation 

4>k^Tk  =  aik  (ti  -  T<u)  6rk. 

This  is  satisfied  by  the  boundary  condition 

4>k  =  o.ik  ( n  -  Tdi)  ■  (46) 

Assuming  arbitrary  variations  in  STk,  this  condition  is  also  necessary. 

In  order  to  control  the  surface  pressure  and  normal  stress  one  can  measure  the  difference 

Tlj  {<Tkj  +  Skj  {p  ~  Pd)}  , 

where  pd  is  the  desired  pressure.  The  normal  component  is  then 

Tn  =  nknj(Tkj  +P~Pd , 

so  that  the  measure  becomes 

*{t)  =  \rl 

=  \ninmnknj  {<xim  +  6lm(p-  pd)}  {crkj  +6kj  (p~  Pd)}  ■ 

JU 

This  corresponds  to  setting 

aik  =  nink 

in  equation  (46).  Defining  the  viscous  normal  stress  as 

Tvn  =  ‘R'k'R'j&kji 

the  measure  can  be  expanded  as 

A f{r)  =  \ninmnknjaim(Tkj  +  ^  (nknj<Tkj  +  nirimcrim)  (p  -  Pd)  +  -  (p  -  Pdf 
=  \Tvn  +  Tvn  (P  -  Pd)  +2  (P-Pdf- 


For  cancellation  of  the  boundary  terms 

<pk  ( njSakj  +  nk6p)  =  {n/nma/m  +  nf  (p  -  Vd)}  nk  ( nj6crkj  +  nk6p) 

leading  to  the  boundary  condition 

<j>k  =  Tlk  (rvn  "bP  Pd)  ■ 

In  the  case  of  high  Reynolds  number,  this  is  well  approximated  by  the  equations 

<t>k  =  nk{jp-  Pd) ,  (47) 

which  should  be  compared  with  the  single  scalar  equation  derived  for  the  inviscid  boundary  condition  (29). 
In  the  ca=r  of  an  inviscid  flow,  choosing 

•^(r)  =  ^  ip  -  Pd)2 

requires 

4>knk5p  =(jp-pd)  nkSp  =  (p-pd)  Sp 

which  is  satisfied  by  equation  (47),  but  which  represents  an  overspecification  of  the  boundary  condition  since 
only  the  single  condition  (29)  need  be  specified  to  ensure  cancellation. 


Boundary  Conditions  Arising  from  the  Energy  Equation 

The  form  of  the  boundary  terms  arising  from  the  energy  equation  depends  on  the  choice  of  temperature 
boundary  condition  at  the  wall.  For  the  adiabatic  case,  the  boundary  contribution  is  (41) 


f  ,  ^dd 

JBk5Td^dB(’ 


while  for  the  constant  temperature  case  the  boundary  term  is  (43).  One  possibility  is  to  introduce  a  contri¬ 
bution  into  the  cost  function  which  depends  on  T  or  so  that  the  appropriate  cancellation  would  occur. 
Since  there  is  little  physical  intuition  to  guide  the  choice  of  such  a  cost  function  for  aerodynamic  design,  a 
more  natural  solution  is  to  set 


0  =  0 


in  the  constant  temperature  case  or 

in  the  adiabatic  case.  Note  that  in  the  constant  temperature  case,  this  choice  of  0  on  the  boundary  would 
also  elimirate  the  boundary  metric  variation  terms  in  (42). 


8  Implementation  of  Navier-Stokes  Design 

The  design  procedures  can  be  summarized  as  follows: 

1.  Solve  the  flow  equations  for  p,  u\,  112,  U3,  P- 

2.  Solve  the  adjoint  equations  for  ip  subject  to  appropriate  boundary  conditions. 

3.  Evaluate  Q  . 

4.  Project  Q  into  an  allowable  subspace  that  satisfies  any  geometric  constraints. 

5.  Update  the  shape  based  on  the  direction  of  steepest  descent. 

6.  Return  to  1  until  convergence  is  reached. 

Practical  implementation  of  the  viscous  design  method  relies  heavily  upon  fast  and  accurate  solvers  for  both 
the  state  (w)  and  co-state  (ip)  systems.  This  work  uses  well-validated  software  for  the  solution  of  the  Euler 
and  Navier-Stokes  equations  developed  over  the  course  of  many  years  [24-26]. 

For  inverse  design  the  lift  is  fixed  by  the  target  pressure.  In  drag  minimization  it  is  also  appropriate  to  fix 
the  lift  coefficient,  because  the  induced  drag  is  a  major  fraction  of  the  total  drag,  and  this  could  be  reduced 
simply  by  reducing  the  lift.  Therefore  the  angle  of  attack  is  adjusted  during  the  flow  solution  to  force  a 
specified  lift  coefficient  to  be  attained,  and  the  influence  of  variations  of  the  angle  of  attack  is  included  in  the 
calculation  of  the  gradient.  The  vortex  drag  also  depends  on  the  span  loading,  which  may  be  constrained  by 
other  considerations  such  as  structural  loading  or  buffet  onset.  Consequently,  the  option  is  provided  to  force 
the  span  loading  by  adjusting  the  twist  distribution  as  well  as  the  angle  of  attack  during  the  flow  solution. 


Discretization 

Both  the  flow  and  the  adjoint  equations  are  discretized  using  a  semi-discrete  cell-centered  finite  volume 
scheme  The  convective  fluxes  across  cell  interfaces  are  represented  by  simple  arithmetic  averages  of  the 
fluxes  computed  using  values  from  the  cells  on  either  side  of  the  face,  augmented  by  artificial  diffusive  terms 
to  prevent  numerical  oscillations  in  the  vicinity  of  shock  waves.  Continuing  to  use  the  summation  convention 
for  repeated  indices,  the  numerical  convective  flux  across  the  interface  between  cells  A  and  B  m  a  three 
dimensional  mesh  has  the  form  ^ 

hAB  =  2^ABj  ttA>  +  ~  ^AB' 

where  SAB  is  the  component  of  the  face  area  in  the  jth  Cartesian  coordinate  direction,  (fAj)  and  (fBj) 
denote  the'flux  fj  as  defined  by  equation  (12)  and  dAB  is  the  diffusive  term.  Variations  of  the  computer 
program  provide  options  for  alternate  constructions  of  the  diffusive  flux. 

The  simplest  option  implements  the  Jameson-Schmidt-Turkel  scheme  [24,27],  using  scalar  diffusive  terms 

of  the  form 

dAB  =  e(2)  Aw  -  e(4)  (Aw+  -  2 Aw  +  Aw  ) , 


where 

Aw  =  wB  —  WA 

and  Aw+  and  Aw~  are  the  same  differences  across  the  adjacent  cell  interfaces  behind  cell  A  and  beyond 
cell  B  in  the  AB  direction.  By  making  the  coefficient  e<2)  depend  on  a  switch  proportional  to  the  undivided 
second  difference  of  a  flow  quantity  such  as  the  pressure  or  entropy,  the  diffusive  flux  becomes  a  third  order 
quantity,  proportional  to  the  cube  of  the  mesh  width  in  regions  where  the  solution  is  smooth.  Oscillations  are 
suppressed  near  a  shock  wave  because  e(2)  becomes  of  order  unity,  while  e(4)  is  reduced  to  zero  by  the  same 
switch.  For  a  scalar  conservation  law,  it  is  shown  in  reference  [27]  that  e(2)  and  e  can  be  constructed  to 
make  the  scheme  satisfy  the  local  extremum  diminishing  (LED)  principle  that  local  maxima  cannot  increase 

while  local  minima  cannot  decrease.  . 

The  second  option  applies  the  same  construction  to  local  characteristic  variables.  There  are  derived  from 

the  eigenvectors  of  the  Jacobian  matrix  Aab  which  exactly  satisfies  the  relation 

Aab  (wb  -  wA)  =  SABj  (fBj  -  fAj)  ■ 


This  corresponds  to  the  definition  of  Roe  [28].  The  resulting  scheme  is  LED  in  the  characteristic  variables. 
The  third  option  implements  the  H-CUSP  scheme  proposed  by  Jameson  [29]  which  combines  differences 
fB  -  fA  and  wB-wA  in  a  manner  suca  that  stationary  shock  waves  can  be  captured  with  a  single  interior 
point  in  the  discrete  solution.  This  scheme  minimizes  the  numerical  diffusion  as  the  velocity  approaches  zero 
in  the  boundary  layer,  and  has  therefore  been  preferred  for  viscous  calculations  in  this  work. 

Similar  artificial  diffusive  terms  are  introduced  in  the  discretization  of  the  adjoint  equation,  but  with 
the  opposite  sign  because  the  wave  directions  are  reversed  in  the  adjoint  equation.  Satisfactory  results  have 
been  obtained  using  scalar  diffusion  in  the  adjoint  equation  while  characteristic  or  H-CUSP  constructions 
are  used  in  the  flow  solution. 

The  discretization  of  the  viscous  terms  of  the  Navier  Stokes  equations  requires  the  evaluation  of  the 
velocity  derivatives  ^  in  order  to  calculate  the  viscous  stress  tensor  Oij  defined  in  equation  (9).  These 

are  most  conveniently  evaluated  at  the  cell  vertices  of  the  primary  mesh  by  introducing  a  dual  mesh  which 
connects  the  cell  centers  of  the  primary  mesh,  as  depicted  in  Figure  (1).  According  to  the  Gauss  formula  for 
a  control  volume  V  with  boundary  S 

[  dv  =  f  UiTljdS 
Jy  OXj  Js 

where  rtj  is  the  outward  normal.  Applied  to  the  dual  cells  this  yields  the  estimate 


dvi 

dxj 


£  6in>s 

faces 


where  S  is  the  area  of  a  face,  and  u*  is  an  estimate  of  the  average  of  over  that  face.  In  order  to  determine 
the  viscous  flux  balance  of  each  primary  cell,  the  viscous  flux  across  each  of  its  faces  is  then  calculated  from 


Fig.  1.  Cell-centered  scheme,  a y  evaluated  at  vertices  of  the  primary  mesh 


the  average  of  the  viscous  stress  tensor  at  the  four  vertices  connected  by  that  face.  This  leads  to  a  compact 
scheme  with  a  stencil  connecting  each  cell  to  its  26  nearest  neighbors. 

The  semi-discrete  schemes  for  both  the  flow  and  the  adjoint  equations  are  both  advanced  to  steady  state 
by  a  multi-stage  time  stepping  scheme.  This  is  a  generalized  Runge-Kutta  scheme  in  which  the  convective  and 
diffusive  terms  are  treated  differently  to  enlarge  the  stability  region  [27,30].  Convergence  to  a  steady  state  is 
accelerated  by  residual  averaging  and  a  multi-grid  procedure  [31].  These  algorithms  have  been  implemented 
both  for  single  and  multiblock  meshes  and  for  operation  on  parallel  computers  with  message  passing  using 
the  MPI  (Message  Passing  Interface)  protocol  [9,32,33]. 

In  this  work,  the  adjoint  and  flow  equations  are  discretized  separately.  The  alternative  approach  of 
deriving  the  discrete  adjoint  equations  directly  from  the  discrete  flow  equations  yields  another  possible 
discretization  of  the  adjoint  partial  differential  equation  which  is  more  complex.  If  the  resulting  equations 
were  solved  exactly,  they  could  provide  the  exact  gradient  of  the  inexact  cost  function  which  results  from  the 
discretization  of  the  flow  equations.  On  the  other  hand,  any  consistent  discretization  of  the  adjoint  partial 
differential  equation  will  yield  the  exact  gradient  as  the  ricsh  is  refined,  and  separate  discretization  has 
proved  to  work  perfectly  well  in  practice.  It  should  also  be  noted  that-  the  discrete  gradient  includes  both 
mesh  effects  and  numerical  errors  such  as  spurious  entropy  production  which  may  not  reflect  the  true  cost 
function  of  the  continuous  problem. 


Mesh  Generation  and  Geometry  Control 

Meshes  for  both  viscous  optimization  and  for  the  treatment  of  complex  configurations  are  externally  gen¬ 
erated  in  order  to  allow  for  their  inspection  and  careful  quality  control.  Single  block  meshes  with  a  C-H 
topology  have  been  used  for  viscous  optimization  of  wing-body  combinations,  while  multiblock  meshes  have 
been  generated  for  complex  configurations  using  GRIDGEN  [34].  In  either  case  geometry  modifications  are 
accommodated  by  a  grid  perturbation  scheme.  For  viscous  wing-body  design  using  single  block  meshes,  the 
wing  surface  mesh  points  themselves  are  taken  as  the  design  variables.  A  simple  mesh  perturbation  scheme 
is  then  used,  in  which  the  mesh  points  lying  on  a  mesh  line  projecting  out  from  the  wing  surface  are  all 
shifted  in  the  same  sense  as  the  surface  mesh  point,  with  a  decay  factor  proportional  to  the  arc  length  along 
the  mesh  line.  The  resulting  perturbation  in  the  face  areas  of  the  neighboring  cells  are  then  included  in 
the  gradient  calculation.  For  complex  configurations  the  geometry  is  controlled  by  superposition  of  analytic 
“bump”  functions  defined  over  the  surfaces  which  are  to  be  modified.  The  grid  is  then  perturbed  to  conform 
to  modifications  of  the  surface  shape  by  the  WARP3D  and  WARP-MB  algorithms  described  in  [32]. 


Optimization 

Two  main  search  procedures  have  been  used  in  our  applications  to  date.  The  first  is  a  simple  descent  method 
in  which  small  steps  are  taken  in  the  negative  gradient  direction.  Let'iF  represent  the  design  variable,  and  Q 
the  gradient.  Then  the  iteration 

87  =  -xg 

can  be  regarded  as  simulating  the  time  dependent  process 


d7 

dt 


=  -G 


where  A  is  the  time  step  At.  Let  A  be  the  Hessian  matrix  with  elements 

_  egt  a2/ 

ij  dTj  dTidTj ' 

Suppose  that  a  locally  minimum  value  of  the  cost  function  I*  =  1(7* )  is  attained  when  7  =  7* .  Then  the 
gradient  G*  =  G(7*)  must  be  zero,  while  the  Hessian  matrix  A *  =  A(7*)  must  be  positive  definite.  Since 
g •  is  zero,  the  cost  function  can  be  expanded  as  a  Taylor  series  in  the  neighborhood  of  7*  with  the  form 

=  r  +  JF*)  +  ... 


Correspondingly, 

G(7)  =  A(7-7*)  +  ... 

As  7  approaches  7* ,  the  leading  terms  become  dominant.  Then,  setting  7  =  (7-7*),  the  search  process 
approximates 

f  =  -A*7. 
dt 

Also,  since  A*  is  positive  definite  it  can  be  expanded  as 

A*  =  RMRt, 

where  M  is  a  diagonal  matrix  containing  the  eigenvalues  of  A *,  and 

RRt  =  RtR  =  L 


Setting 


v  =  Rt  T, 


the  search  process  can  be  represented  as 


dv 

dt 


—Mv. 


The  stability  region  for  the  simple  forward  Euler  stepping  scheme  is  a  unit  circle  centered  at  -1  on  the 
negative  real  axis.  Thus  for  stability  we  must  choose 


M max^  ~  A*max^  <  2, 


while  the  asymptotic  decay  rate,  given  by  the  smallest  eigenvalue,  is  proportional  to 

In  order  to  make  sure  that  each  new  shape  in  the  optimization  sequence  remains  smooth,  it  proves 
essential  to  smooth  the  gradient  and  to  replace  G  by  its  smoothed  value  G  in  the  descent  process.  This  also 
acts  as  a  preconditioner  which  allows  the  use  of  much  larger  steps.  To  apply  smoothing  in  the  £1  direction, 
for  example,  the  smoothed  gradient  G  may  be  calculated  from  a  discrete  approximation  to 


(48) 


where  e  is  the  smoothing  parameter.  If  one  sets  ST  =  -A Q,  then,  assuming  the  modification  is  applied  on 
the  surface  £,2  =  constant,  the  first  order  change  in  the  cost  function  is 

SI  =  -  J J  QSTd£,\d£,i 

=  -A//(«2  +  .(g))«A 
<0, 

assuring  an  improvement  if  A  is  sufficiently  small  and  positive,  unless  the  process  has  already  reached  a 
stationary  point  at  which  £5  =  0. 

It  turns  out  that  this  approach  is  tolerant  to  the  use  of  approximate  values  of  the  gradient,  so  that 
neither  the  flow  solution  nor  the  adjoint  solution  need  be  fully  converged  before  making  a  shape  change. 
This  results  in  very  large  savings  in  the  computational  cost.  For  inviscid  optimization  it  is  necessary  to  use 
only  15  multigrid  cycles  for  the  flow  solution  and  the  adjoint  solution  in  each  design  iteration.  For  viscous 

optimization,  about  20-30  multigrid  cycles  are  needed. 

Our  second  main  search  procedure  incorporates  a  quasi-Newton  method  for  general  constrained  opti¬ 
mization.  In  this  class  of  methods  the  step  is  defined  as 

ST  =  -XPQ, 

where  P  is  a  preconditioner  for  the  search.  An  ideal  choice  is  P  =  A*  1,  so  that  the  corresponding  time 
dependent  process  reduces  to 


for  which  all  the  eigenvalues  are  equal  to  unity,  and  f  is  reduced  to  zero  in  one  time  step  by  the  choice 
At  =  1  if  the  Hessian,  A,  is  constant.  The  full  Newton  method  takes  P  =  A"1,  requiring  the  evaluation  of  the 
Hessian  matrix,  A,  at  each  step.  It  corresponds  to  the  use  of  the  Newton-Raphson  method  to  solve  the  non¬ 
linear  equation  Q  =  0.  Quasi-Newton  methods  estimate  A*  from  the  change  in  the  gradient  during  the  search 
process.  This  requires  accurate  estimates  of  the  gradient  at  each  time  step.  In  order  to  obtain  these,  both 
the  flow  solution  and  the  adjoint  equation  must  be  fully  converged.  Most  quasi-Newton  methods  also  require 
a  line  search  in  each  search  direction,  for  which  the  flow  equations  and  cost  function  must  be  accurately 
evaluated  several  times.  They  have  proven  quite  robust  in  practice  for  aerodynamic  optimization  [35]. 

In  the  applications  to  complex  configurations  presented  below  the  optimization  was  carried  out  using 
the  existing,  well  validated  software  NPSOL.  This  software,  which  implements  a  quasi-Newton  method  for 
optimization  with  both  linear  and  non-linear  constraints,  has  proved  very  reliable  but  is  generally  more 
expensive  than  the  simple  search  method  with  smoothing. 


9  Industrial  Experience  and  Results 

The  methods  described  in  this  paper  have  been  quite  thoroughly  tested  in  industrial  applications  in  which 
they  were  used  as  a  tool  for  aerodynamic  design.  They  have  proved  useful  both  in  inverse  mode  to  find 
shapes  that  would  produce  desired  pressure  distributions,  and  for  direct  minimization  of  the  drag.  They 
have  been  applied  both  to  well  understood  configurations  that  have  gradually  evolved  through  incremental 
improvements  guided  by  wind  tunnel  tests  and  computational  simulation,  and  to  new  concepts  for  which 
there  is  a  limited  knowledge  base.  In  either  case  they  have  enabled  engineers  to  produce  improved  designs. 

Substantial  improvements  are  usually  obtained  with  20  —  200  design  cycles,  depending  on  the  difficulty  of 
the  case.  One  concern  is  the  possibility  of  getting  trapped  in  a  local  minimum.  In  practice  this  has  not  proved 
to  be  a  source  of  difficulty.  In  inverse  mode,  it  often  proves  possible  to  come  very  close  to  realizing  the  target 
pressure  distribution,  thus  effectively  demonstrating  convergence.  In  drag  minimization,  the  result  of  the 
optimization  is  usually  a  shock-free  wing.  If  one  considers  drag  minimization  of  airfoils  in  two-dimensional 


inviscid  transonic  flow,  it  can  be  seen  that  every  shock-free  airfoil  produces  zero  drag,  and  thus  optimization 
based  solely  on  drag  has  a  highly  non-unique  solution.  Different  shock-free  airfoils  can  be  obtained  by  starting 
from  different  initial  profiles.  One  may  also  influence  the  character  of  the  final  design  by  blending  a  target 
pressure  distribution  with  the  drag  in  the  definition  of  the  cost  function. 

Similar  considerations  apply  to  three-dimensional  wing  design.  Since  the  vortex  drag  can  be  reduced 
simply  by  reducing  the  lift,  the  lift  coefficient  must  be  fixed  for  a  meaningful  drag  minimization.  In  order 
to  do  this  the  angle  of  attack  a  is  adjusted  during  the  flow  solution.  It  has  proved  most  effective  to  make 
a  small  change  6a  proportional  to  the  difference  between  the  actual  and  the  desired  lift  coefficient  at  every 
iteration  in  the  flow  calculation.  A  typical  wing  of  a  transport  aircraft  is  designed  for  a  lift  coefficient  in 
the  range  of  0.4  to  0.6.  The  total  wing  drag  may  be  broken  down  into  vortex  drag,  drag  due  to  viscous 
effects,  and  shock  drag.  The  vortex  drag  coefficient  is  typically  in  the  range  of  0.0100  (100  counts),  while 
the  friction  drag  coefficient  is  in  the  range  of  45  counts,  and  the  shock  drag  at  a  Mach  number  just  before 
the  onset  of  severe  drag  rise  is  of  the  order  of  15  counts.  With  a  fixed  span,  typically  dictated  by  structural 
limits  or  a  constraint  imposed  by  airport  gates,  the  vortex  drag  is  entirely  a  function  of  span  loading,  and  is 
minimized  by  an  elliptic  loading  unless  winglets  are  added.  Transport  aircraft  usually  have  highly  tapered 
wings  with  very  large  root  chords  to  accommodate  retraction  of  the  undercarriage.  An  elliptic  loading  may 
lead  to  excessively  large  section  lift  coefficients  on  the  outboard  wing,  leading  to  premature  shock  stall  or 
buffet  when  the  load  is  increased.  The  structure  weight  is  also  reduced  by  a  more  inboard  loading  which 
reduces  the  root  bending  moment.  Thus  the  choice  of  span  loading  is  influenced  by  other  considerations.  The 
skin  friction  of  transport  aircraft  is  typically  very  close  to  flat  plate  skin  friction  in  turbulent  flow,  and  is  very 
insensitive  to  section  variations.  An  exception  to  this  is  the  case  of  smaller  executive  jet  aircraft,  for  which 
the  Reynolds  number  may  be  small  enough  to  allow  a  significant  run  of  laminar  flow  if  the  suction  peak 
of  the  pressure  distribution  is  moved  back  on  the  section.  This  leaves  the  shock  drag  as  the  primary  target 
for  wing  section  optimization.  This  is  reduced  to  zero  if  the  wing  is  shock-free,  leaving  no  room  for  further 
improvement.  Thus  the  attainment  of  a  shock-free  flow  is  a  demonstration  of  a  successful  drag  minimization. 
In  practice  range  is  maximized  by  maximizing  M  ,  and  this  is  likely  to  be  increased  by  increasing  the  lift 
coefficient  to  the  point  where  a  weak  shock  appears.  One  may  also  use  optimization  to  find  the  maximum 
Mach  number  at  which  the  shock  drag  can  be  eliminated  or  significantly  reduced  for  a  wing  with  a  given 
sweepback  angle  and  thickness.  Alternatively  one  may  try  to  find  the  largest  wing  thickness  or  the  minimum 
sweepback  angle  for  which  the  shock  drag  can  be  eliminated  at  a  given  Mach  number.  This  can  yield  both 
savings  in  structure  weight  and  increased  fuel  volume  .  If  there  is  no  fixed  limit  for  the  wing  span,  such  as 
a  gate  constraint,  increased  thickness  can  be  used  to  allow  an  increase  in  aspect  ratio  for  a  wing  of  equal 
weight,  in  turn  leading  to  a  reduction  in  vortex  drag.  Since  the  vortex  drag  is  usually  the  largest  component 
of  the  tofai  wing  drag,  this  is  probably  the  most  effective  design  strategy,  and  it  may  pay  to  increase  the 
wing  thickness  to  the  point  where  the  optimized  section  produces  a  weak  shock  wave  rather  than  a  shock-free 
flow  [23]. 

The  first  major  industrial  application  of  an  adjoint  based  aerodynamic  optimization  method  was  the 
wing  design  of  the  Beech  Premier  [36]  in  1995.  The  method  was  successfully  used  in  inverse  mode  as  a  tool 
to  obtain  pressure  distributions  favorable  to  the  maintenance  of  natural  laminar  flow  over  a  range  of  cruise 
Mach  numbers.  Wing  contours  were  obtained  which  yielded  the  desired  pressure  distribution  in  the  presence 
of  closely  coupled  engine  nacelles  on  the  fuselage  above  the  wing  trailing  edge. 

During  1996  some  preliminary  studies  indicated  that  the  wings  of  both  the  McDonnell  Douglas  MD-11 
and  the  Boeing  747-200  could  be  made  shock-free  in  a  representative  cruise  condition  by  using  very  small 
shape  modifications,  with  consequent  drag  savings  which  could  amount  to  several  percent  of  the  total  drag. 
This  led  to  a  decision  to  evaluate  adjoint-based  design  methods  in  the  design  of  the  McDonnell  Douglas 
MDXX  during  the  summer  and  fall  of  1996.  In  initial  studies  wing  redesigns  were  carried  out  for  inviscid 
transonic  flow  modeled  by  the  Euler  equations.  A  redesign  to  minimize  the  drag  at  a  specified  lift  and  Mach 
number  required  about  40  design  cycles,  which  could  be  completed  overnight  on  a  workstation. 

Three  main  lessons  were  drawn  from  these  initial  studies:  (i)  the  fuselage  effect  is  to  large  to  be  ignored 
and  must  be  included  in  the  optimization,  (ii)  single-point  designs  could  be  too  sensitive  to  small  variations 
in  the  flight  condition,  typically  producing  a  shock-free  flow  at  the  design  point  with  a  tendency  to  break  up 
into  a  severe  double  shock  pattern  below  the  design  point,  and  (iii)  the  shape  changes  necessary  to  optimize 
a  wing  in  transonic  flow  are  smaller  than  the  boundary  layer  displacement  thickness,  with  the  consequence 
that  viscous  effects  must  be  included  in  the  final  design. 


In  order  to  meet  the  first  two  of  these  considerations,  the  second  phase  of  the  study  was  concentrated 
on  the  optimization  of  wing-body  combinations  with  multiple  design  points.  These  were  still  performed  with 
inviscid  flow  to  reduce  computational  cost  and  allow  for  fast  turnaround.  It  was  found  that  comparatively 
insensitive  designs  could  be  obtained  by  minimizing  the  drag  at  a  fixed  Mach  number  for  three  fairly  closely 
spaced  lift  coefficients  such  as  0.5,  0.525,  and  0.55,  or  alternatively  three  nearby  Mach  numbers  with  a  fixed 
lift  coefficient. 

The  third  phase  of  the  project  was  focused  on  the  design  with  viscous  effects  using  as  a  starting  point 
wings  which  resulted  from  multipoint  inviscid  optimization.  While  the  full  viscous  adjoint  method  was  still 
under  development,  it  was  found  that  useful  improvements  could  be  realized,  particularly  in  inverse  mode, 
using  the  inviscid  result  to  provide  the  target  pressure,  by  coupling  an  inviscid  adjoint  solver  to  a  viscous 
flow  solver.  Computer  costs  are  many  times  larger,  both  because  finer  meshes  are  needed  to  resolve  the 
boundary  layer,  and  because  more  iterations  are  needed  in  the  flow  and  adjoint  solutions.  In  order  to  force 
the  specified  lift  coefficient  the  number  of  iterations  in  each  flow  solution  had  to  be  increased  from  15  to 
100.  To  achieve  overnight  turnaround  a  fully  parallel  implementation  of  the  software  had  to  be  developed. 
Finally  it  was  found  that  in  order  to  produce  sufficiently  accurate  results,  the  number  of  mesh  points  had  to 
be  increased  to  about  1.8  million.  In  the  final  phase  of  this  project  it  was  planned  to  carry  out  a  propulsion 
integration  study  using  the  multiblock  versions  of  the  software.  This  study  was  not  completed  due  to  the 
cancellation  of  the  entire  MDXX  project. 

In  the  next  subsections  we  present  examples  of  the  use  of  the  adjoint  method  for  viscous  inverse  and 
drag  minimization  in  two  dimensional  flow.  We  then  show  a  three-dimensional  wing  design  using  the  Euler 
equations  and  a  wing  design  using  the  full  viscous  adjoint  method  in  its  current  form,  implemented  in  the 
computer  program  SYN107.-  These  calculations  were  all  performed  using  the  simple  descent  method  with 
smoothing  of  the  gradient.  This  has  proved  to  be  very  efficient:  in  all  cases  the  final  optimum  design  was 
achieved  with  a  total  computational  cost  equivalent  to  the  cost  of  from  2  to  10  converged  flow  solutions.  The 
remaining  subsections  present  results  of  optimizations  for  complete  configurations  in  inviscid  transonic  and 
supersonic  flow  using  the  multiblock  parallel  design  program,  SYN107-MB. 

Inverse  design  of  an  airfoil  in  transonic  viscous  flow 

Our  first  example  shows  an  inverse  design  in  two  dimensional  viscous  transonic  flow  obtained  using  the 
two-dimensional  design  code  SYN103.  The  target  pressure  is  that  of  the  section  of  the  ONERA  M6  wing 
at  Mach  .75  and  a  lift  coefficient  of  .50.  It  was  calculated  using  SYN103  in  analysis  mode,  thus  it  should 
be  exactly  realizable.  A  C-type  mesh  was  used  which  contained  256  intervals  in  the  chordwise  and  96  cells 
in  the  normal  direction  for  a  total  of  24, 576  cells.  The  design  calculation  was  started  with  the  NACA  0012 
airfoil  as  the  initial  profile,  and  tbs  ONERA  M6  pressure  distribution  was  almost  exactly  recovered  in  25 
design  cycles.  In  the  first  cycle  120  iterations  were  used  in  both  the  flow  and  the  adjoint  solutions.  In  the 
subsequent  cycles  only  30  iterations  were  used  in  both  the  flow  and  adjoint  solutions.  Figure  2  shows 
the  initial  profile  and  pressure  distribution  with  the  pressure  coefficient  plotted  vertically  in  the  negative 
direction.  It  then  shows  the  results  after  one,  five  and  twenty  five  design  cycles,  with  the  target  represented 
by  circles.  It  also  superposes  on  each  redesigned  profile  the  smoothed  gradient  plotted  in  the  direction  of  the 
shape  modification.  A  fixed  scale  is  used  so  that  it  is  possible  to  observe  the  decrease  in  the  magnitude  of 
the  gradient  as  the  calculation  converges  enough  to  ensure  that  they  were  fairly  close  to  convergence.  The 
root  mean  square  error  between  the  target  and  actual  pressure  was  reduced  from  .0530  to  .0016  in  the  course 
of  the  entire  calculation  which  took  3569  seconds  using  a  single  Silicon  Graphics  R10000  processor.  A  fully 
converged  flow  solution  using  500  iterations  on  the  same  mesh  took  936  seconds,  so  the  cost  of  the  entire 
design  calculation  was  about  that  of  three  flow  solutions. 


Drag  reduction  of  an  airfoil  in  viscous  flow 

The  next  example  shows  a  redesign  of  the  RAE2822  airfoil  to  reduce  the  drag  at  a  fixed  lift  coefficient  of 
.65  in  transonic  flow  at  Mach  .75.  In  this  case  a  shock  free  flow  was  obtained  after  10  design  cycles,  in 
each  of  which  both  the  flow  and  the  adjoint  solutions  were  calculated  with  25  multigrid  cycles.  A  grid  with 
512  x  64  cells  was  used.  The  pressure  drag  was  reduced  from  .0091  to  .0041,  while  the  viscous  drag  remained 
essentially  constant.  The  constraint  was  imposed  that  the  thickness  of  the  profile  could  not  be  reduced 


bv  only  permitting  outward  movement  from  the  initial  profile.  Figure  3  displays  the  sequence  of  pressure 
distributions  showing  the  elimination  of  the  shock  wave.  It  also  shows  the  initial  profile,  and  the  smoothed 
gradient  superposed  In  the  subsequent  profiles.  It  can  be  «„  that  the  gradient  cont:  nues »  have  an inward 
component  indicating  that  the  drag  might  be  further  reduced  if  a  thickness  reduction  were  permitted  It 
should  be  noted  that  the  unsmoothed  gradient  is  in  the  sense  of  crossing  over  the  trailing  edge,  because  the 
resulting  non  physical  shape  would  correspond  to  a  sink  in  the  free  stream  which  would  have  a  negative 
drag.  The  solution  of  the  smoothing  equation  (48)  with  a  two  point  boundary  condition  allow  the  trailing 

edge  to  be  frozen. 

Three  point  inviscid  redesign  of  the  Boeing  747  wing 

The  third  example  shows  a  redesign  of  the  wing  of  the  Boeing  747  to  reduce  its  drag  in  a  typical  cruising 
condition.  It  has  been  our  experience  that  drag  minimization  at  a  single  pomt  tends  to  produce  a  wing 
which  is  shock  free  at  its  design  point,  but  tends  to  display  undesirable  characteristics  off  its  design  point. 
Typically,  a  double  shock  pattern  forms  below  the  design  lift  coefficient  and  Mach  number,  and  a  single 
fairly  strong  shock  above  the  design  point.  To  alleviate  this  tendency  the  calculation  was  perfonned  wi 
three  design  points.  In  carrying  out  multipoint  designs  of  this  kind  a  composite  gradient  is  calculated  as  a 
weighted  average  of  the  gradients  calculated  for  each  design  point  separately.  In  this  case  the  design  points 
were  selected  as  lift  coefficients  of  .38, .42  and  .46  for  the  exposed  wing  at  Mach  .85.  Because  the  fuselage  has 
a  significant  effect  on  the  flow  over  the  wing,  the  calculations  were  performed  for  the  wing  body  combmation 
buAhe  shape  modifications  were  restricted  to  the  wing  alone.  The  fuselage  also  contributes  to  the  lift,  so 
that  the  total  lift  coefficient  at  the  mid  design  point  was  estimated  to  be  .50. 

The  results  are  displayed  in  Figures  4-  6  and  in  Table  1  which  shows  the  drag  at  three  design  points 
of  the  initial  wing,  and  the  final  wing  after  30  design  cycles.  It  can  been  seen  that  a  drag  reduction  was 
obtained  over  the  entire  range  of  lift  coefficients,  and  at  the  mid  design  point  the  redesigned  wing  is  almost 
shock  free  Figure  7  shows  the  modification  in  the  wing  section  about  half  way  out  the  span.  It  can  be  seen 
that  a  useful  drag  reduction  can  be  obtained  by  a  very  small  change  in  the  wing  shape.  Thls  ls  bec^e  ° 
the  extreme  sensitivity  of  the  transonic  flow.  Also,  it  is  clear  that  without  a  tool  of  this  kind  it  would  be 

almost  impossible  to  find  an  optimum  shape. 


Design  Conditions  Initial  Three  Point  Design 

Mach  Cl  Cd  Original  Cp  Redesign  Cp  Reduction  (%) 


0.85 

0.38 

0.0071 

0.0064 

9.8 

0.85 

0.42 

0.0086 

0.0U7 

10.4 

0.85 

0.46 

0.0106 

0  OOdS 

10.3 

Table  1.  Drag  Reduction  for  Multipoint  Design. 


Transonic  Viscous  Wing-Body  Design 

A  typical  result  for  drag  minimization  of  a  wing  body  combination  in  transonic  viscous  flow  is  presented 
next  The  viscous  adjoint  optimization  method  was  used  with  a  Baldwin-Lomax  turbulence  model.  The 
initial  wing  is  similar  to  one  produced  during  the  MDXX  design  studies.  Figures  8-10  show  the  result  of  the 
wing-body  redesign  on  a  C-H  mesh  with  288  x  96  x  64  cells.  The  wing  has  sweep  back  of  about  38  degrees 
at  the  1/4  chord.  A  total  of  44  iterations  of  the  viscous  optimization  procedure  resulted  m  a  shock-free  wing 
at  a  cruise  design  point  of  Mach  0.86,  with  a  lift  coefficient  of  0.61  for  the  wing-body  combmation  at  a 
Reynolds  number  of  101  million  based  on  the  root  chord.  Using  48  processors  of  an  SGI  Ongm2000  parallel 
computer,  each  design  iteration  takes  about  22  minutes  so  that  overnight  turnaround  for  such  a  calculation 
is  possible.  Figure  8  compares  the  pressure  distribution  of  the  final  design  with  that  of  the  initial  wmg.  e 
final  wing  is  quite  thick,  with  a  thickness  to  chord  ratio  of  about  14  percent  at  the  root  and  9  percent  at 
the  tip.  The  optimization  was  performed  with  a  constraint  that  the  section  modifications  were  not  allowe 


to  decrease  the  thickness  anywhere.  The  design  offers  excellent  performance  at  the  nominal  cruise  point.  A 
drag  reduction  of  2.2  counts  was  achieved  from  the  initial  wing  which  had  itself  been  derived  by  inviscid 
optimization.  Figures  9  and  10  show  the  results  of  a  Mach  number  sweep  to  determine  the  dragnse.  The 
drag  coefficients  shown  in  the  figures  represent  the  total  wing  drag  including  shock,  vortex,  and  skin  friction 
contributions.  It  can  be  seen  that  a  double  shock  pattern  forms  below  the  design  point,  while  there  is  actually 
a  slight  increase  in  the  drag  coefficient  at  Mach  0.85.  This  wing  has  a  low  drag  coefficient  however,  over  a 
wide  range  of  conditions.  Above  the  design  point  a  single  shock  forms  and  strengthens  as  the  Mach  number 

increases. 


10  Conclusions 

The  adjoint  design  method  developed  under  this  grant  is  now  well  established  and  has  proved  effective  in  a 
variety  of  industrial  applications.  The  method  combines  the  versatility  of  numerical  optimization  methods 
with  the  efficiency  of  inverse  design  methods.  The  geometry  is  modified  by  a  pid  perturbation  technique 
which  is  applicable  to  arbitrary  configurations.  Both  the  wing-body  and  multiblock  version  of  the  design 
algorithms  have  been  implemented  in  parallel  using  the  MPI  (Message  Passing  Interface)  Standard,  and 
they  both  yield  excellent  parallel  speedups.  The  combination  of  computational  efficiency  with  geometric 
flexibility  provides  a  powerful  tool,  with  the  final  goal  being  to  create  practical  aerodynamic  shape  design 
methods  for  complete  aircraft  configurations. 
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2a:  Cp  after  Zero  Design  Cycles. 
Design  Mach  0.75t  Ci  =  0.5008, 
Cd  =  0.0225. 
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2b:  Cp  after  One  Design  Cycle. 
Design  Mach  0.75,  Ci  =  0.4841, 
Cd  =  0.0185. 
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2c:  Cp  after  Five  Design  Cycles. 
Design  Mach  0.75,  Cj  =  0.4994, 
Cd  =  0.0148. 


2d:  Cp  after  Twenty  five  Design  Cycles. 
Design  Mach  0.75,  C/  =  0.5007, 

Cd  =  0.0118. 


Fig.  2.  Inverse  Design  of  an  ONERA  Airfoil.  The  vectors  on  the  airfoil  surface  represent  the  direction  and  magnitude 
of  the  gradient. 


3a:  Cp  after  Zero  Design  Cycles. 
Design  Mach  0.75,  C\  =  0.6450, 
Cd(pressure)  =  0.0091,' 
Cd(viscous)  =  0.0056. 
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3b:  Cp  after  One  Design  Cycle. 
Design  Mach  0.75,  C<  =  0.6512, 
Cd(pressure)  =  0.0066. 
Cd(m5Coii5)  —  0.0057. 


3c:  Cp  after  Two  Design  Cycles. 
Design  Mach  0.75,  Cj  =  0.6510, 
C'd  (pressure)  =  0.0054, 

Cd  (viscous)  =  0.0057. 


3d:  Cp  after  Ten  Design  Cycles. 
Design  Mach  0.75,  Cj  =  0.6460, 
Cd{pressure)  —  0.0041. 
Cd(viscous)  =  0.0058. 


Fig.  3.  Drag  Minimization  of  an  RAE2822  Airfoil.  The  vectors  on  the  airfoil  surface  represent  the  direction  and 
magnitude  of  the  gradient. 
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Fig.  4.  Pressure  distribution  of  the  Boeing  747  Wing-Body  before  optimization. 


COMPARISON  OF  CHORDWISE  PRESSURE  DISTRIBUTIONS 
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Fig.  5.  Pressure  distribution  of  the  Boeing  747  Wing-Body  after  a  three  point  optimization. 
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COMPARISON  OF  CHORDWISE  PRESSURE  DISTRIBUTIONS 
MPX5X  WING-BODY  \* 


Fig.  8.  Pressure  distribution  of  the  MPX5X  before  and  after  optimization. 
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COMPARISON  OF  CHORDWISE  PRESSURE  DISTRIBUTIONS 
MPX5X  WING-BODY 


Fig.  10.  Off  design  performance  of  the  MPX5X  above  the  design  point. 


